Data access and discovery

Authors & Contributors

Authors

• Pier Lorenzo Marasco, Ispra (Italy), @pl-marasco

• Alejandro Coca-Castro, The Alan Turing Institute (United Kingdom), @acocac

Contributors

• Tina Odaka, Ifremer (France), @tinaok

Overview

Questions

• How to access online (remote) datasets?
• How to discover online geoscience data?
• What is openeo?
• What is stac?
• Why do chunking and compressing matter?
• What is Analysis Ready, cloud optimized data (ARCO)?
• What is pangeo-forge?

Objectives

• Learn to access datasets from online object storage
• Learn about online catalogs
• Learn about openeo
• Learn about stac
• Learn about chunking and compression
• Learn about Analysis Cloud Optimized data (ARCO)
• Learnn about Pangeo Forge initiative

Context

We will be using Long Term Statistics (1999-2019) product provided by the Copernicus Global Land Service over Lombardia and access them through S3-comptabile storage (OpenStack Object Storage “Swift”) and openeo Sentinel-2 Cloud-Optimised Dataset online through STAC.

Data

This episode we will use data available online or that has been prepared and stored in Zenodo.

Setup

This episode uses the following Python packages:

• xarray

• netcdf4

• h5netcdf

• wget

• numpy

• s3fs

Please install these packages if not already available in your Python environment. Below, we only install packages that are not available in the EGI-ACE EOSC deployment of Pangeo for the FOSS4G course.

Packages

In this episode, Python packages are imported when we start to use them. However, for best software practices, we recommend you to install and import all the necessary libraries at the top of your Jupyter notebook.

Introduction to the Long Term statistics

CGLS LTS are computed over a time span of 20 years aggregated over each 10 days period (month/01,month/11, month/21). For each date the long term minimum, maximum, mean, median and standard deviation are computed.

S3-compatible Object Storage to access online data

Up to now we downloaded data locally and then open the downloaded file with Xarray open_dataset. Each attendee downloaded a local copy locally. When willing to manipulate large amount of data, this approach is not optimal (requires a lot of unecessary local downloads). Sharing data online as Object Storage allows for data sharing and access to much larger amount of data.

One of the most popular way to access online remote data is through Amazon Simple Storage Service (S3) and you do not necessarily need to use Amazon service to benefit from S3 Object storage. Many other providers offers S3-compatible object storage that can be access in a very similar way.

Below we will be accessing online the NDVI Long Term statistics from Copernicus Land Service that we have publicly stored in OpenStack Object storage (swift).

import s3fs
import xarray as xr

fs = s3fs.S3FileSystem(anon=True,
      client_kwargs={
         'endpoint_url': 'https://object-store.cloud.muni.cz'
      })

Tip

The parameter anon is for anonymous and is set to True because the data we have stored at https://object-store.cloud.muni.cz is public

List files and folders in existing buckets

Instead of organizing files in various folders, object storage systems store files in a flat organization of containers (called “buckets”).

fs.ls('foss4g-data')

['foss4g-data/CGLS_LTS_1999_2019',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia',
 'foss4g-data/test']

fs.ls('foss4g-data/CGLS_LTS_1999_2019_Lombardia')

['foss4g-data/CGLS_LTS_1999_2019_Lombardia/',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0101_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0111_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0121_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0201_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0211_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0221_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0301_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0311_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0321_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0401_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0411_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0421_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0501_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0511_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0521_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0601_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0611_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0621_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0701_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0711_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0721_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0801_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0811_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0821_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0901_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0911_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0921_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1001_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1011_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1021_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1101_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1111_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1121_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1201_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1211_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1221_GLOBE_VGT-PROBAV_V3.0.1.nc']

Access remote files from S3-compatible Object Storage

s3path = 's3://foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0101_GLOBE_VGT-PROBAV_V3.0.1.nc'

LTS = xr.open_dataset(fs.open(s3path))

LTS

<xarray.Dataset>
Dimensions:  (lon: 331, lat: 235)
Coordinates:
  * lon      (lon) float64 8.5 8.509 8.518 8.527 ... 11.42 11.43 11.44 11.45
  * lat      (lat) float64 46.7 46.69 46.68 46.67 ... 44.63 44.63 44.62 44.61
Data variables:
    crs      |S1 b''
    min      (lat, lon) float32 nan nan -0.076 -0.048 ... 0.216 0.256 0.236
    median   (lat, lon) float32 nan nan -0.072 -0.048 ... 0.404 0.404 0.4 0.416
    max      (lat, lon) float32 nan nan -0.028 -0.048 ... 0.56 0.54 0.568 0.588
    mean     (lat, lon) float32 nan nan -0.064 -0.048 ... 0.388 0.392 0.404 0.42
    stdev    (lat, lon) float32 nan nan 0.024 nan ... 0.088 0.088 0.088 0.104
    nobs     (lat, lon) float32 nan nan 4.0 1.0 nan ... 21.0 21.0 21.0 21.0 21.0
Attributes: (12/19)
    Conventions:          CF-1.6
    parent_identifier:    urn:cgls:global:ndvi_stats_all
    identifier:           urn:cgls:global:ndvi_stats_all:NDVI-LTS_1999-2019-0...
    long_name:            Normalized Difference Vegetation Index
    title:                Normalized Difference Vegetation Index: Long Term S...
    product_version:      V3.0.1
    ...                   ...
    source:               Derived from EO satellite imagery
    processing_mode:      Offline
    references:           https://land.copernicus.eu/global/products/ndvi
    copyright:            Copernicus Service information 2021
    archive_facility:     VITO
    history:              2021-03-01 - Processing line NDVI LTS

xarray.Dataset

• Dimensions:
  • lon: 331
  • lat: 235
• Coordinates: (2)
  • lon
    (lon)
    float64
    8.5 8.509 8.518 ... 11.44 11.45

    _CoordinateAxisType :

    Lon

    axis :

    X

    DIMENSION_LABELS :

    lon

    long_name :

    longitude

    standard_name :

    longitude

    units :

    degrees_east

    array([ 8.5     ,  8.508929,  8.517857, ..., 11.428571, 11.4375  , 11.446429])

  • lat
    (lat)
    float64
    46.7 46.69 46.68 ... 44.62 44.61

    _CoordinateAxisType :

    Lat

    axis :

    Y

    DIMENSION_LABELS :

    lat

    long_name :

    latitude

    standard_name :

    latitude

    units :

    degrees_north

    array([46.696429, 46.6875  , 46.678571, ..., 44.625   , 44.616071, 44.607143])

• Data variables: (7)
  • crs
    ()
    |S1
    ...

    long_name :

    coordinate reference system

    _CoordinateAxisTypes :

    GeoX GeoY

    _CoordinateTransformType :

    Projection

    grid_mapping_name :

    latitude_longitude

    inverse_flattening :

    298.257223563

    longitude_of_prime_meridian :

    0.0

    semi_major_axis :

    6378137.0

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -180.0000000000 0.0089285714 0.0 80.0000000000 0.0 -0.0089285714

    array(b'', dtype='|S1')

  • min
    (lat, lon)
    float32
    ...

    standard_name :

    normalized_difference_vegetation_index

    long_name :

    Minimum Normalized Difference Vegetation Index over time period

    grid_mapping :

    crs

    flag_meanings :

    sea no_data

    flag_values :

    [254 255]

    units :

    valid_range :

    [ 0 250]

    cell_methods :

    area: mean time: minimum

    array([[   nan,    nan, -0.076, ...,  0.216,  0.192,  0.208],
           [   nan,    nan,    nan, ...,  0.196,  0.204,  0.244],
           [ 0.124,    nan,    nan, ...,  0.068,  0.048,  0.048],
           ...,
           [ 0.244,  0.24 ,  0.224, ...,  0.176,  0.22 ,  0.232],
           [ 0.144,  0.128,  0.232, ...,  0.216,  0.256,  0.24 ],
           [ 0.18 ,  0.18 ,  0.16 , ...,  0.216,  0.256,  0.236]], dtype=float32)

  • median
    (lat, lon)
    float32
    ...

    standard_name :

    normalized_difference_vegetation_index

    long_name :

    Median Normalized Difference Vegetation Index over time period

    grid_mapping :

    crs

    flag_meanings :

    sea no_data

    flag_values :

    [254 255]

    units :

    valid_range :

    [ 0 250]

    cell_methods :

    area: mean time: median

    array([[   nan,    nan, -0.072, ...,  0.26 ,  0.26 ,  0.26 ],
           [   nan,    nan,    nan, ...,  0.216,  0.296,  0.332],
           [ 0.124,    nan,    nan, ...,  0.188,  0.264,  0.264],
           ...,
           [ 0.392,  0.384,  0.388, ...,  0.372,  0.392,  0.392],
           [ 0.352,  0.372,  0.376, ...,  0.392,  0.412,  0.424],
           [ 0.4  ,  0.404,  0.424, ...,  0.404,  0.4  ,  0.416]], dtype=float32)

  • max
    (lat, lon)
    float32
    ...

    standard_name :

    normalized_difference_vegetation_index

    long_name :

    Maximum Normalized Difference Vegetation Index over time period

    grid_mapping :

    crs

    flag_meanings :

    sea no_data

    flag_values :

    [254 255]

    units :

    valid_range :

    [ 0 250]

    cell_methods :

    area: mean time: maximum

    array([[   nan,    nan, -0.028, ...,  0.332,  0.3  ,  0.344],
           [   nan,    nan,    nan, ...,  0.272,  0.324,  0.404],
           [ 0.124,    nan,    nan, ...,  0.288,  0.464,  0.44 ],
           ...,
           [ 0.52 ,  0.52 ,  0.484, ...,  0.496,  0.532,  0.528],
           [ 0.488,  0.48 ,  0.488, ...,  0.544,  0.568,  0.56 ],
           [ 0.5  ,  0.592,  0.52 , ...,  0.54 ,  0.568,  0.588]], dtype=float32)

  • mean
    (lat, lon)
    float32
    ...

    standard_name :

    normalized_difference_vegetation_index

    long_name :

    Mean Normalized Difference Vegetation Index over time period

    grid_mapping :

    crs

    flag_meanings :

    sea no_data

    flag_values :

    [254 255]

    units :

    valid_range :

    [ 0 250]

    cell_methods :

    area: mean time: mean

    array([[   nan,    nan, -0.064, ...,  0.268,  0.252,  0.272],
           [   nan,    nan,    nan, ...,  0.224,  0.28 ,  0.328],
           [ 0.124,    nan,    nan, ...,  0.2  ,  0.276,  0.272],
           ...,
           [ 0.384,  0.38 ,  0.376, ...,  0.344,  0.38 ,  0.384],
           [ 0.34 ,  0.352,  0.364, ...,  0.376,  0.416,  0.42 ],
           [ 0.384,  0.396,  0.388, ...,  0.392,  0.404,  0.42 ]], dtype=float32)

  • stdev
    (lat, lon)
    float32
    ...

    standard_name :

    normalized_difference_vegetation_index

    long_name :

    Standard deviation Normalized Difference Vegetation Index over time period

    grid_mapping :

    crs

    flag_meanings :

    sea no_data

    flag_values :

    [254 255]

    units :

    valid_range :

    [ 0 250]

    cell_methods :

    area: mean time: standard_deviation

    array([[  nan,   nan, 0.024, ..., 0.048, 0.052, 0.068],
           [  nan,   nan,   nan, ..., 0.036, 0.052, 0.08 ],
           [  nan,   nan,   nan, ..., 0.068, 0.144, 0.152],
           ...,
           [0.08 , 0.08 , 0.072, ..., 0.088, 0.096, 0.096],
           [0.092, 0.096, 0.072, ..., 0.096, 0.1  , 0.096],
           [0.088, 0.104, 0.104, ..., 0.088, 0.088, 0.104]], dtype=float32)

  • nobs
    (lat, lon)
    float32
    ...

    standard_name :

    normalized_difference_vegetation_index number_of_observations

    grid_mapping :

    crs

    units :

    valid_range :

    [ 1 250]

    cell_methods :

    area: mean

    array([[nan, nan,  4., ...,  4.,  4.,  3.],
           [nan, nan, nan, ...,  4.,  4.,  4.],
           [ 1., nan, nan, ...,  9.,  6.,  5.],
           ...,
           [19., 20., 20., ..., 21., 21., 21.],
           [20., 18., 19., ..., 20., 21., 21.],
           [20., 19., 18., ..., 21., 21., 21.]], dtype=float32)

• Attributes: (19)

  Conventions :

  CF-1.6

  parent_identifier :

  urn:cgls:global:ndvi_stats_all

  identifier :

  urn:cgls:global:ndvi_stats_all:NDVI-LTS_1999-2019-0101_GLOBE_V3.0.1

  long_name :

  Normalized Difference Vegetation Index

  title :

  Normalized Difference Vegetation Index: Long Term Statistics 1KM: GLOBE 1999-2019 0101

  product_version :

  V3.0.1

  time_coverage_start :

  1999-01-01T00:00:00Z

  time_coverage_end :

  2019-12-31T23:59:59Z

  platform :

  SPOT-4, SPOT-5, Proba-V

  sensor :

  VEGETATION-1, VEGETATION-2, VEGETATION

  orbit_type :

  LEO

  processing_level :

  L4

  institution :

  VITO NV

  source :

  Derived from EO satellite imagery

  processing_mode :

  Offline

  references :

  https://land.copernicus.eu/global/products/ndvi

  copyright :

  Copernicus Service information 2021

  archive_facility :

  VITO

  history :

  2021-03-01 - Processing line NDVI LTS

LTS.sel(lat=45.88, lon=8.63, method='nearest')['min'].values

array(0.264, dtype=float32)

Warning

The same dataset can be available from different locations e.g. CGLS distributor VITO, Zenodo, S3-compatible OpenStack Object storage (Swift), etc. How do you know if the analyze the very same dataset? You cannot know except if the datasets have a persistent identifier such as a Digital Object Identifier. It is therefore recommended to be extra careful with where you get your datasets and double check that the content is what you expect (for instance, you can perform basic quality checks).

Access multiple remote files

s3path = 's3://foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_*.nc'

remote_files = fs.glob(s3path)
remote_files

['foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0101_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0111_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0121_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0201_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0211_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0221_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0301_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0311_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0321_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0401_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0411_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0421_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0501_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0511_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0521_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0601_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0611_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0621_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0701_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0711_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0721_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0801_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0811_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0821_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0901_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0911_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-0921_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1001_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1011_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1021_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1101_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1111_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1121_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1201_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1211_GLOBE_VGT-PROBAV_V3.0.1.nc',
 'foss4g-data/CGLS_LTS_1999_2019_Lombardia/AOI_c_gls_NDVI-LTS_1999-2019-1221_GLOBE_VGT-PROBAV_V3.0.1.nc']

We need to add a time dimension to concatenate data. For this, we define a function that will be called for each remote file (via the preprocess parameter of Xarray open_mfdataset.)

from datetime import datetime

def preprocess(ds):
    t = datetime.strptime(ds.attrs['identifier'].split(':')[-1].split('_')[1].replace('1999-', ''), "%Y-%m%d")
    return(ds.assign_coords(time=t).expand_dims('time'))   

Xarray open_mfdataset allows to open multiple files at the same time.

# Iterate through remote_files to create a fileset
fileset = [fs.open(file) for file in remote_files]

When opening remote files, you can also select the variables you wish to analyze.

LTS = xr.open_mfdataset(fileset, combine='nested', concat_dim=['time'],  preprocess=preprocess,
                        decode_coords="all")
LTS

<xarray.Dataset>
Dimensions:  (lon: 331, lat: 235, time: 36)
Coordinates:
  * lon      (lon) float64 8.5 8.509 8.518 8.527 ... 11.42 11.43 11.44 11.45
  * lat      (lat) float64 46.7 46.69 46.68 46.67 ... 44.63 44.63 44.62 44.61
    crs      |S1 b''
  * time     (time) datetime64[ns] 2019-01-01 2019-01-11 ... 2019-12-21
Data variables:
    min      (time, lat, lon) float32 dask.array<chunksize=(1, 235, 331), meta=np.ndarray>
    median   (time, lat, lon) float32 dask.array<chunksize=(1, 235, 331), meta=np.ndarray>
    max      (time, lat, lon) float32 dask.array<chunksize=(1, 235, 331), meta=np.ndarray>
    mean     (time, lat, lon) float32 dask.array<chunksize=(1, 235, 331), meta=np.ndarray>
    stdev    (time, lat, lon) float32 dask.array<chunksize=(1, 235, 331), meta=np.ndarray>
    nobs     (time, lat, lon) float32 dask.array<chunksize=(1, 235, 331), meta=np.ndarray>
Attributes: (12/19)
    Conventions:          CF-1.6
    parent_identifier:    urn:cgls:global:ndvi_stats_all
    identifier:           urn:cgls:global:ndvi_stats_all:NDVI-LTS_1999-2019-0...
    long_name:            Normalized Difference Vegetation Index
    title:                Normalized Difference Vegetation Index: Long Term S...
    product_version:      V3.0.1
    ...                   ...
    source:               Derived from EO satellite imagery
    processing_mode:      Offline
    references:           https://land.copernicus.eu/global/products/ndvi
    copyright:            Copernicus Service information 2021
    archive_facility:     VITO
    history:              2021-03-01 - Processing line NDVI LTS

xarray.Dataset

• Dimensions:
  • lon: 331
  • lat: 235
  • time: 36
• Coordinates: (4)
  • lon
    (lon)
    float64
    8.5 8.509 8.518 ... 11.44 11.45

    _CoordinateAxisType :

    Lon

    axis :

    X

    DIMENSION_LABELS :

    lon

    long_name :

    longitude

    standard_name :

    longitude

    units :

    degrees_east

    array([ 8.5     ,  8.508929,  8.517857, ..., 11.428571, 11.4375  , 11.446429])

  • lat
    (lat)
    float64
    46.7 46.69 46.68 ... 44.62 44.61

    _CoordinateAxisType :

    Lat

    axis :

    Y

    DIMENSION_LABELS :

    lat

    long_name :

    latitude

    standard_name :

    latitude

    units :

    degrees_north

    array([46.696429, 46.6875  , 46.678571, ..., 44.625   , 44.616071, 44.607143])

  • crs
    ()
    |S1
    b''

    long_name :

    coordinate reference system

    _CoordinateAxisTypes :

    GeoX GeoY

    _CoordinateTransformType :

    Projection

    grid_mapping_name :

    latitude_longitude

    inverse_flattening :

    298.257223563

    longitude_of_prime_meridian :

    0.0

    semi_major_axis :

    6378137.0

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -180.0000000000 0.0089285714 0.0 80.0000000000 0.0 -0.0089285714

    array(b'', dtype='|S1')

  • time
    (time)
    datetime64[ns]
    2019-01-01 ... 2019-12-21

    array(['2019-01-01T00:00:00.000000000', '2019-01-11T00:00:00.000000000',
           '2019-01-21T00:00:00.000000000', '2019-02-01T00:00:00.000000000',
           '2019-02-11T00:00:00.000000000', '2019-02-21T00:00:00.000000000',
           '2019-03-01T00:00:00.000000000', '2019-03-11T00:00:00.000000000',
           '2019-03-21T00:00:00.000000000', '2019-04-01T00:00:00.000000000',
           '2019-04-11T00:00:00.000000000', '2019-04-21T00:00:00.000000000',
           '2019-05-01T00:00:00.000000000', '2019-05-11T00:00:00.000000000',
           '2019-05-21T00:00:00.000000000', '2019-06-01T00:00:00.000000000',
           '2019-06-11T00:00:00.000000000', '2019-06-21T00:00:00.000000000',
           '2019-07-01T00:00:00.000000000', '2019-07-11T00:00:00.000000000',
           '2019-07-21T00:00:00.000000000', '2019-08-01T00:00:00.000000000',
           '2019-08-11T00:00:00.000000000', '2019-08-21T00:00:00.000000000',
           '2019-09-01T00:00:00.000000000', '2019-09-11T00:00:00.000000000',
           '2019-09-21T00:00:00.000000000', '2019-10-01T00:00:00.000000000',
           '2019-10-11T00:00:00.000000000', '2019-10-21T00:00:00.000000000',
           '2019-11-01T00:00:00.000000000', '2019-11-11T00:00:00.000000000',
           '2019-11-21T00:00:00.000000000', '2019-12-01T00:00:00.000000000',
           '2019-12-11T00:00:00.000000000', '2019-12-21T00:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (6)
  • min
    (time, lat, lon)
    float32
    dask.array<chunksize=(1, 235, 331), meta=np.ndarray>

    standard_name :

    normalized_difference_vegetation_index

    long_name :

    Minimum Normalized Difference Vegetation Index over time period

    flag_meanings :

    sea no_data

    flag_values :

    [254 255]

    units :

    valid_range :

    [ 0 250]

    cell_methods :

    area: mean time: minimum

    Array Chunk Bytes 10.68 MiB 303.85 kiB Shape (36, 235, 331) (1, 235, 331) Count 109 Graph Layers 36 Chunks Type float32 numpy.ndarray

  • median
    (time, lat, lon)
    float32
    dask.array<chunksize=(1, 235, 331), meta=np.ndarray>

    standard_name :

    normalized_difference_vegetation_index

    long_name :

    Median Normalized Difference Vegetation Index over time period

    flag_meanings :

    sea no_data

    flag_values :

    [254 255]

    units :

    valid_range :

    [ 0 250]

    cell_methods :

    area: mean time: median

    Array Chunk Bytes 10.68 MiB 303.85 kiB Shape (36, 235, 331) (1, 235, 331) Count 109 Graph Layers 36 Chunks Type float32 numpy.ndarray

  • max
    (time, lat, lon)
    float32
    dask.array<chunksize=(1, 235, 331), meta=np.ndarray>

    standard_name :

    normalized_difference_vegetation_index

    long_name :

    Maximum Normalized Difference Vegetation Index over time period

    flag_meanings :

    sea no_data

    flag_values :

    [254 255]

    units :

    valid_range :

    [ 0 250]

    cell_methods :

    area: mean time: maximum

    Array Chunk Bytes 10.68 MiB 303.85 kiB Shape (36, 235, 331) (1, 235, 331) Count 109 Graph Layers 36 Chunks Type float32 numpy.ndarray

  • mean
    (time, lat, lon)
    float32
    dask.array<chunksize=(1, 235, 331), meta=np.ndarray>

    standard_name :

    normalized_difference_vegetation_index

    long_name :

    Mean Normalized Difference Vegetation Index over time period

    flag_meanings :

    sea no_data

    flag_values :

    [254 255]

    units :

    valid_range :

    [ 0 250]

    cell_methods :

    area: mean time: mean

    Array Chunk Bytes 10.68 MiB 303.85 kiB Shape (36, 235, 331) (1, 235, 331) Count 109 Graph Layers 36 Chunks Type float32 numpy.ndarray

  • stdev
    (time, lat, lon)
    float32
    dask.array<chunksize=(1, 235, 331), meta=np.ndarray>

    standard_name :

    normalized_difference_vegetation_index

    long_name :

    Standard deviation Normalized Difference Vegetation Index over time period

    flag_meanings :

    sea no_data

    flag_values :

    [254 255]

    units :

    valid_range :

    [ 0 250]

    cell_methods :

    area: mean time: standard_deviation

    Array Chunk Bytes 10.68 MiB 303.85 kiB Shape (36, 235, 331) (1, 235, 331) Count 109 Graph Layers 36 Chunks Type float32 numpy.ndarray

  • nobs
    (time, lat, lon)
    float32
    dask.array<chunksize=(1, 235, 331), meta=np.ndarray>

    standard_name :

    normalized_difference_vegetation_index number_of_observations

    units :

    valid_range :

    [ 1 250]

    cell_methods :

    area: mean

    Array Chunk Bytes 10.68 MiB 303.85 kiB Shape (36, 235, 331) (1, 235, 331) Count 109 Graph Layers 36 Chunks Type float32 numpy.ndarray

• Attributes: (19)

  Conventions :

  CF-1.6

  parent_identifier :

  urn:cgls:global:ndvi_stats_all

  identifier :

  urn:cgls:global:ndvi_stats_all:NDVI-LTS_1999-2019-0101_GLOBE_V3.0.1

  long_name :

  Normalized Difference Vegetation Index

  title :

  Normalized Difference Vegetation Index: Long Term Statistics 1KM: GLOBE 1999-2019 0101

  product_version :

  V3.0.1

  time_coverage_start :

  1999-01-01T00:00:00Z

  time_coverage_end :

  2019-12-31T23:59:59Z

  platform :

  SPOT-4, SPOT-5, Proba-V

  sensor :

  VEGETATION-1, VEGETATION-2, VEGETATION

  orbit_type :

  LEO

  processing_level :

  L4

  institution :

  VITO NV

  source :

  Derived from EO satellite imagery

  processing_mode :

  Offline

  references :

  https://land.copernicus.eu/global/products/ndvi

  copyright :

  Copernicus Service information 2021

  archive_facility :

  VITO

  history :

  2021-03-01 - Processing line NDVI LTS

Tip

If you use one of xarray’s open methods such as xarray.open_dataset to load netCDF files with the default engine, it is recommended to use decode_coords=”all”. This will load the grid mapping variable into coordinates for compatibility with rioxarray. See rioxarray documentation.

Chunking and compression

You may have missed it but if you look carefully to LTS, each Data Variable is now a dask.array with a chunsize of (1, 15680, 40320). So basically accessing one data variable would load (1, 15680, 40320) into the memory (and you need quite a large amount of memory to load such dataset!).

For netCDF files, the chunks correspond to the entire size of the dataset. We will see later that other data formats such as zarr are much more optimized for accessing and analyzing large amount of data.

Select Area Of Interest

Clip data according to a polygon

We use the same shapefile as the one used in the previous episode.

import geopandas as gpd

try:
    GAUL = gpd.read_file('Italy.geojson')
except:
    GAUL = gpd.read_file('zip+https://mars.jrc.ec.europa.eu/asap/files/gaul1_asap.zip') 

AOI_name = 'Lombardia'
AOI = GAUL[GAUL.name1 == AOI_name]
AOI_poly = AOI.geometry
AOI_poly

14    POLYGON ((10.23973 46.62177, 10.25084 46.61110...
Name: geometry, dtype: geometry

Define the reference system

LTS = LTS.rio.write_crs(4326)

Once this has been done we can clip data with the polygon that has been read through geopandas at the beginning of the notebook.

LTS = LTS.rio.clip(AOI_poly, crs=4326)

Save Clipped LTS into local zarr

LTS.to_zarr(rf'c_gls_NDVI-LTS_1999-2019_AOI_VGT-PROBAV_V3.zarr', mode='w')

<xarray.backends.zarr.ZarrStore at 0x7f58520830d0>

Read zarr from local file

LTS = xr.open_zarr(rf'c_gls_NDVI-LTS_1999-2019_AOI_VGT-PROBAV_V3.zarr')

Visualize LTS with Holoviews

# visualization
import matplotlib.pyplot as plt
import holoviews as hv
import hvplot.xarray
import hvplot.pandas

LTS['min'][5].hvplot(geo=True, cmap="RdYlGn")

LTS_NDVI_POI = LTS.sel(lat=45.88, lon=8.63, method='nearest')

POI_min = LTS_NDVI_POI['min']
POI_max = LTS_NDVI_POI['max']
POI_std = LTS_NDVI_POI['stdev']
POI_mean = LTS_NDVI_POI['mean']

import numpy as np

(
 POI_mean.hvplot(c='grey', label='mean LTS 1999-2019', alpha=0.2, line_dash='dashed') *
 POI_min.hvplot(c='b',label='min LTS 1999-2019', alpha=0.2, line_dash='dashed') *
 POI_max.hvplot(c='r', label='Max LTS 1999-2019', alpha=0.2, line_dash='dashed') *
 hv.Area((POI_mean.time, POI_mean - 2 * POI_std, POI_mean + 2 * POI_std ), vdims=['- Std', '+ Std'], label='Standard deviation LTS 1999-2019',).opts( color='blue', alpha=0.2))\
    .opts(title=f"CGLS S3 300m NDVI fluctuation over the year 2022 in comparison with CGLS Long Term Statistics\nPoint of interest Lat {np.round(POI.lat.data,2)} Lon {np.round(POI.lon.data,2)}",
          legend_position='bottom_right')

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [27], in <cell line: 2>()
      1 (
      2  POI_mean.hvplot(c='grey', label='mean LTS 1999-2019', alpha=0.2, line_dash='dashed') *
      3  POI_min.hvplot(c='b',label='min LTS 1999-2019', alpha=0.2, line_dash='dashed') *
      4  POI_max.hvplot(c='r', label='Max LTS 1999-2019', alpha=0.2, line_dash='dashed') *
      5  hv.Area((POI_mean.time, POI_mean - 2 * POI_std, POI_mean + 2 * POI_std ), vdims=['- Std', '+ Std'], label='Standard deviation LTS 1999-2019',).opts( color='blue', alpha=0.2))\
----> 6     .opts(title=f"CGLS S3 300m NDVI fluctuation over the year 2022 in comparison with CGLS Long Term Statistics\nPoint of interest Lat {np.round(POI.lat.data,2)} Lon {np.round(POI.lon.data,2)}",
      7           legend_position='bottom_right')

NameError: name 'POI' is not defined

Discover online datasets

There are so many available online datasets that it can be very difficult to find out where to find the data you need. To ease your work, data catalogs have been defined and allow to browse and search datasets

Spatio Temporal Asset Catalogs (STAC)

The STAC specification is a common language to describe geospatial information, so it can more easily be worked with, indexed, and discovered.

Openeo

openEO develops an open API to connect R, Python, JavaScript and other clients to big Earth observation cloud back-ends in a simple and unified way.

Analysis Ready, cloud optimized data (ARCO)

Accessing LTS from OpenStack Object Storage through a catalog

The Pangeo forge initiative

Key Points

• Online catalogs
• STAC
• Openeo
• Chunking and compression
• Pangeo Forge